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Abstract 

The  relativistic  radiation  damping  and  accompanying  electron  self-force  is,  in  principle,  properly 
modeled  by  the  standard  PIC  (particle-in-cell)  algorithm;  however,  due  to  grid  resolution  and 
time  step  constraints,  there  is  no  possibility  of  directly  simulating  the  fields  of  the  emitted 
radiation  over  a  spectrum  sufficient  to  include  a  significant  fraction  of  the  radiated  power,  so  a 
scheme  has  been  developed  to  correct  for  these  errors.  A  model  for  secondary  emission  of 
electrons  from  arbitrary  surfaces  with  user  definable  coefficients  has  been  developed  and 
implemented  in  the  XOOPIC  (XI 1 -based  Object-Oriented  Particle-In-Cell)  code,  allowing  for 
separate  secondary  emission  treatment  of  each  species  incident  on  each  boundary.  A  literature 
survey  determined  the  four  algorithms  most  relevant  to  treating  electromagnetic  boundary 
conditions  in  a  PIC  model,  along  with  their  respective  strengths  and  weaknesses.  The  errors  in 
the  generic  PIC  algorithm  have  been  analyzed  in  the  electrostatic  limit,  and  an  initial  extension  of 
this  analysis  to  the  generic  electromagnetic  PIC  scheme  was  developed.  A  large  subset  of  the 
existing  XOOPIC  GUI,  which  runs  only  under  Unix,  has  been  ported  to  a  new  cross-platform 
windowing  library,  called  Qt,  that  runs  equally  well  on  Microsoft  Windows  95/98/NT  and  all 
versions  of  Unix. 


Modeling  the  Electron  Self-Force 

The  relativistic  electron  self  force  is  a  physical  process  in  which  the  acceleration  of  a  relativistic 
electron  causes  the  electron  to  radiate.  The  energy  is  negligible  except  for  radiation  from  highly 
relativistic  particles.  A  prototypical  configuration  that  exhibits  radiation  damping  is  a  relativistic 
electron  in  a  strong  magnetic  field. 

Jackson  1  derives  the  power  radiated  from  a  relativistic  particle  undergoing  a  strong  acceleration 
in  a  uniform  magnetic  field: 
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where  e  is  the  electronic  charge,  c  is  the  speed  of  light,  e0  is  the  permittivity  of  free  space, 
P  =  v/c ,  and  y  =  ( 1  —  >82 ) 1/2 .  The  relativistic  gyroradius  is  given  by: 
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where  m  is  the  electron  mass  and  B  is  the  magnetic  field. 
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The  energy  radiated  by  the  electron  per  gyro-orbit  is  approximately: 

p  _  ^  «3„,4 
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The  spectral  content  of  the  radiated  power  is  a  complicated  function.  Jackson  derives  the  energy 
radiated  per  unit  frequency  per  unit  solid  angle: 
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Here,  6  is  the  angle  from  the  plane  of  circular  motion,  K  is  the  modified  Bessel  function,  and  t, 
is  given  by: 
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The  radiation  is  strongly  polarized  in  the  plane  of  motion,  with  increasing  confinement  to  the 
plane  of  motion  with  increasing  frequency.  The  radiation  is  negligible  beyond  the  critical 
frequency. 
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Another  useful  quantity  is  the  critical  harmonic  number,  ncril  =  3 y3 ,  which  describes  the  number 

of  significant  harmonics  present  in  the  radiation  spectrum  for  a  particle  undergoing  a  periodic 
motion  such  as  synchrotron  radiation. 


Similarly,  a  critical  angle  beyond  which  the  radiation  drops  to  1/e  is  given  by: 
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Equation  6  indicates  that  the  radiated  power  is  confined  to  increasingly  smaller  angles  at  higher 
frequencies,  as  well  as  with  increasing  particle  energy.  The  radiation  is  emitted  from  a  particle 
undergoing  circular  motion  in  a  pattern  similar  to  the  sweep  of  a  lighthouse  beam,  with  the  focus 
always  centered  in  the  direction  of  the  velocity. 

The  relativistic  radiation  damping  and  accompanying  self  force  is  properly  modeled  by  the 
coupled  Maxwell  and  Newton-Lorentz  Equations.  However,  in  PIC  simulation,  the  finite  nature 
of  the  time-  and  space-steps  results  in  error  for  this  force.  An  example  is  useful  in  illustrating  the 
difficulty  in  modeling  this  phenomena  with  the  standard  PIC  method. 

Consider  a  1  GeV  electron  orbiting  a  10  T  static  magnetic  field.  The  relativistic  momentum 
factor  is  7  =  1952,  and  the  velocity  is  effectively  the  speed  of  light,  c,  to  six  places.  The 
gyroradius  is  0.334  m.  The  power  radiated  is  P  =  6xl0'6  W,  so  the  electron  would  lose  1%  of 
its  energy  in  about  270  ns.  The  critical  frequency  is  cocril  I 2k  =  3.2xl018  Hz,  with  a 
corresponding  critical  wavelength  of  Xcril  =9.4xl0“u  m.  In  order  to  resolve  the  wavelength  on 
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the  mesh,  we  would  need  Ax  <  A/4  =  2.4x1 0"11  m.  For  a  nominal  system  length  of  L-  3 rL  in 
each  dimension,  the  required  number  of  cells  in  each  dimension  is  4.2xl010,  or  1.7xl021  cells 
in  a  two  dimensional  simulation.  This  exceeds  the  memory  capacity  of  existing  computers. 
Furthermore,  the  timestep  required  to  satisfy  the  Courant  condition  is  At  <  Ax/V2 c  =  5.7  xlO-20 
s.  The  cyclotron  period  is  7  ns,  so  one  orbit  would  require  1.2x10"  timesteps,  while  running  to 
1%  energy  loss  would  require  4.7 xlO12  timesteps.  Resolving  the  critical  frequency  requires  a 
timestep  At  <  3.1xl(T19  s,  which  is  even  more  untenable. 

The  conclusion  to  this  analysis  is  that  the  radiation  damping  cannot  be  self-consistently  modeled 
by  the  standard  PIC  scheme  even  for  a  single  particle.  Due  to  the  grid  resolution  and  timestep 
constraints,  there  is  no  possibility  of  directly  simulating  the  fields  of  the  radiation  emitted  over  a 
spectrum  sufficient  to  include  a  significant  fraction  of  the  radiated  power. 

To  correct  this  problem,  one  can  implement  a  modification  to  the  equation  of  motion  that  will 
include  the  radiation  damping  term  for  the  particle  equations  of  motion: 

Frad  =  2  e2  7  fdp  dp)  ? 

3  47T£0c4  -Jy2  - 1  \  dt  dt 

where  p  is  the  particle  momentum.  The  derivative  must  be  time  centered  to  maintain  the  second 
order  accuracy  of  the  equation  of  motion: 


dt  At 


Note  that  the  relativistic  momenta  are  computed  at  the  half  timestep  level  in  the  standard  PIC 
leapfrog  mover,  so  this  term  maintains  the  existing  accuracy  of  the  method  (see  accuracy  section 
below).  The  force  is  directed  anti-parallel  to  the  velocity  of  the  electron. 

If  it  is  desired  to  follow  the  radiation  emitted  (for  example  to  measure  the  Poynting  flux  in  some 
solid  angle),  one  might  apply  conservation  of  energy  in  and  invert  the  frequency  and  angular 
distributions  above  to  emit  photons.  The  photons  would  statistically  represent  the  radiation. 
Depending  on  the  spectral  and  angular  resolution  requirements,  this  technique  may  also  result  in 
an  intractable  problem. 

Secondary  Emission  of  Electrons 

A  model  for  secondary  emission  from  arbitrary  surfaces  with  user  definable  coefficients  has  been 
developed  and  implemented  in  the  XOOPIC  code.  The  secondary  coefficient  has  both  angular 
and  energy  dependence,  following  the  model  of  Vaughan2: 

a(y,0)=Sn axo  lH,f^exP(l-w))\  9 

where  energy  of  the  incident  particle  is  denoted  by  V  ,  and  the  angle  of  incidence  measured  from 
the  normal  is  6 .  Here,  <5max0  is  the  peak  in  the  secondary  emission  curve,  which  occurs  for  the 

incident  energy  .  The  parameter,  k,  is  given  by: 
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The  energy  dependence  appears  implicitly  in  Eq.  9,  through  the  normalized  energy  w. 

V-V0 

Here,  V0  is  the  threshold  energy  below  which  the  secondary  coefficient  becomes  zero  (i.e. 
negative  values  of  secondary  coefficient  are  not  allowed). 

The  angular  dependence  appears  explicitly  in  the  equation  for  the  secondary  coefficient  as  well  as 
implicitly  in  the  energy  dependent  term.  The  parameter,  0 < ks  <2,  describes  the  smoothness  of 

the  surface,  with  0  being  rough  and  2  being  smooth.  Vaughan  uses  separate  smoothness  factors  in 
the  angular  and  energy  dependent  parts  of  the  equation,  but  says  that  values  of  1  are 
recommended  unless  data  for  a  specific  surface  supports  other  values. 


Figure  1.  Secondary  electron  energy  emission  distribution 

The  emission  spectrum  of  secondary  electrons  has  three  energy  regimes,  as  shown  in  Figure  1. 
Incident  electrons  reflected  at  the  surface  of  the  secondary  emitting  material  are  called  reflected 
primaries.  Reflected  primaries  comprise  about  3%  of  the  emitted  electron  population.  The  energy 
of  the  reflected  primary  is  approximately  the  same  as  the  energy  of  the  incident  primary,  Vr  =Vr 
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The  primary  is  reflected  specularly  in  angle,  0r  =  -0, ,  where  0  is  measured  from  the  normal  to 
the  surface  of  impact. 

Backscattered  primaries  are  electrons  that  impact  the  surface,  and  scatter  off  of  several  lattice 
atoms  and/or  impurities.  Typically,  backscattered  primaries  comprise  about  7%  of  the  emitted 
electron  population.  These  electrons  re-emerge  from  the  impact  surface  with  energies  in  the 
range  0  <  Vb  <  V. .  Within  this  energy  band,  all  energies  are  taken  as  equally  probable,  Vb  =  RVt , 

where  0  <  R  <  1  is  a  uniformly  distributed  random  number.  The  angle  of  emission  is  specular  just 
as  in  the  reflected  primary  case.  A  more  detailed  treatment  might  consider  a  distribution  of  angles 
resulting  from  quantum  mechanical  treatment  of  the  lattice  potentials  in  the  secondary  emission 
medium. 


True  secondaries  are  electrons  that  are  emitted  from  the  lattice  at  impact  surface.  The  true 
secondaries  comprise  about  90%  of  the  emitted  electron  population.  Energy  is  imparted  to  the 
lattice  electrons  over  some  time  long  compared  to  an  elastic  collision  time,  so  the  electron  energy 
distribution  can  be  approximated  as  a  Maxwell-Boltzmann  distribution.  The  electrons  are  emitted 
with  energies  from  a  thermal  distribution  of  temperature  V, , 


f(y)= 


-v  /  v, 


12 


Due  to  the  timescale  of  the  emission  process,  the  angle  of  emission  is  taken  to  be  isotropic: 

g(0)  =  cos(0)/2.  13 

In  the  XOOPIC  implementation,  all  the  parameters  above  can  be  modified  from  the  input  file, 
with  the  typical  values  above  used  as  defaults.  The  generic  model  can  thus  be  used  to 
approximately  parameterize  secondary  emission  from  many  different  media,  and  indeed  has  been 
used  to  model  electron-induced  secondary  emission  as  well  as  ion-induced  secondary  emission 
common  to  discharges.  In  addition,  the  present  XOOPIC  model  allows  for  separate  secondary 
emission  treatment  for  each  species  incident  on  each  boundary.  For  example,  one  electrode  could 
have  one  profile  for  electrons  and  another  for  ions,  while  another  electrode  might  have  the  same 
secondary  emission  parameters  for  both  species. 


Electromagnetic  Boundary  Conditions 

A  literature  survey  has  been  performed  to  acquire  information  about  electromagnetic  boundary 
conditions.  Here  we  summarize  the  wave  absorbing  methods  most  relevant  to  PIC  modeling, 
including  surface  impedance  boundary  conditions  (SIBC),  the  Lindman  EM  Boundary  Condition, 
Superabsorption  Boundary  Condition,  and  Exact  Non-reflecting  Boundary  Conditions. 

Surface  impedance  boundary  conditions  (SIBC)  (see,  for  example,  Ref.'s  3  and  6)  are  based  on 
the  definition  of  surface  impedances  for  plane  waves  with  horizontal  or  vertical  polarizations. 
This  method  is  employed  to  reduce  the  solution  volume  in  the  analysis  of  the  interactions  of 
electromagnetic  waves  with  objects  of  arbitrary  shape  and  material  composition  (including  free 
space).  Implementation  with  constant  surface  impedance  is  straightforward,  but  it  is  not  accurate 
for  waves  with  a  large  angle  of  incidence.  This  method  can  be  expanded  to  include  dispersive 
surface  impedance. 
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The  Lindman  Boundary  Condition  is  based  primarily  on  the  paper  by  Lindman.4  Propagating  and 
evanescent  waves  from  within  the  computational  region  generates  no  reflected  waves  as  they 
cross  the  Lindman  boundary.  The  boundary  condition  is  formulated  for  a  wave  equation.  An 
alternate  picture  is  to  consider  the  boundary  condition  an  expansion  of  the  normal  impedance 
boundary  condition  to  (almost)  any  angle  of  incidence.  This  boundary  condition  is 
straightforward  to  formulate  with  the  frequency  and  wavenumber  domain,  but  is  difficult  to  use 
in  time  and  space  domain  code.  The  difficulty  lies  in  finding  a  numerically  efficient 
approximation.  A  disadvantage  of  this  method  is  that  there  is  a  transient  response  to  an  incoming 
wave.  The  closer  the  incidence  wave  is  to  90  degree  incidence,  the  longer  the  transient.  Another 
problem  is  that  the  method  is  unstable  when  there  is  a  density  jump  parallel  to  and  near  the 
boundary,  such  as  with  a  plasma  sheath  or  beam  edge. 

The  Superabsorption  Boundary  Layer  can  be  applied  to  known  absorbing  boundary  conditions, 
resulting  in  dramatic  reduction  in  the  numerical  error  caused  by  the  boundary  reflection.  The 
principle  and  analysis  of  the  superabsorption  method  is  presented  in  Mei  and  Fang.5  This  method 
makes  use  of  the  fact  that  when  an  absorbing  boundary  condition  is  applied  to  the  boundary 
electric  fields,  which  are  then  used  to  find  the  magnetic  fields  in  their  neighborhood,  the 
reflection  errors  in  the  magnetic  fields  are  related  to  those  magnetic  fields,  which  are  computed 
directly  from  the  absorbing  boundary  condition.  The  interesting  part  of  the  relation  between  these 
two  errors  is  that  they  bear  the  opposite  sign.  The  results  of  the  two  different  calculations  of  the 
boundary  magnetic  fields  can  be  combined  in  such  a  manner  that  the  errors  cancel,  but  not  the 
fields  remain. 


The  Surface  Integral  Boundary  Condition  is  an  exact  global  method  for  achieving  a  non¬ 
reflecting  boundary  condition  (NRBC).  It  is  based  on  a  time  domain  integral  representation  of  the 
fields  outside  the  calculation  volume  in  terms  of  the  known  fields  on  the  surface  surrounding  that 
volume.  The  robustness  and  accuracy  of  this  method  is  countered  by  the  complexity  and  large 
computational  expense. 

Let  V  be  a  volume  with  a  surface  S.  The  interior  of  V  must  be  homogeneous  and  isotropic:  the 
permittivity,  e ,  and  permeability,  fi ,  are  constant  throughout  the  volume.  Assuming  that  the 
fields  satisfy  the  Sommerfeld  radiation  condition  and  the  source  current  density,  J,  and  charge 
density,  p ,  are  located  in  V,  the  electric  field  at  r  is  given  by: 


\dV 

V 

'  p{r',T)t 

1  i  a,p( t',t)F 

V 

£ 

4 neR1 

— f- 

4k eRc 

4  kR  J 

rfnxafpH(r',T)  fix  1 

(E(r ',  t)  +  3  f  pE(r t  )x  R ) 

]  h-(E(r,,T)+raipE(r,,r)R) 

{  4neR 

s\ 

4neR2 

4 neR2 

/ 

14 


Where  T  =  t-R/c,  the  time  it  takes  a  signal  to  travel  from  r  to  r  at  the  speed  of  light, 
c  =  \/yfeji . 


Analysis  of  Errors  Inherent  in  the  Generic  PIC  Algorithm 

The  errors  generic  PIC  algorithm  have  been  analyzed  separately  in  the  electrostatic  limit  by 
Birdsall  and  Langdon6  and  Hockney  and  Eastwood.7  These  results  can  be  used  to  analyze  the 
error  in  the  entire  PIC  scheme  by  examining  how  errors  can  propagate  from  one  algorithm  to  the 
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next.  To  demonstrate  the  analysis  methodology,  the  electrostatic  PIC  method  is  examined  first. 
Note  that  this  step  is  relevant  in  electromagnetic  schemes  because  the  same  algorithms  are 
employed  to  obtain  the  initial  fields  for  an  initial  configuration  of  charges.  Next,  the  same 
analysis  is  extended  to  the  generic  electromagnetic  PIC  scheme. 

Figure  2  shows  the  schematic  of  error  propagation  in  the  generic  electrostatic  PIC  algorithm. 
Errors  propagate  from  particle  position  to  charge  density  via  interpolation  of  PIC  particle  charges 
to  the  mesh.  The  resulting  error  in  charge  density  is  then  propagated  to  the  Poisson  Equation.  The 
error  in  potential  is  next  propagated  to  the  electric  field  as  a  result  of  the  differencing  scheme 
employed.  The  error  in  the  electric  field  on  the  mesh  is  then  propagated  to  the  Lorentz  Equations 
as  a  result  of  interpolation,  first  to  particle  velocity  and  then  to  the  particle  position.  If  a  static 
magnetic  field  must  also  be  interpolated  from  values  known  on  a  mesh,  then  an  additional  error 
is  introduced  into  the  Lorentz  Equations.  The  error  in  the  Lorentz  Equations  is  then  propagated 
back  to  the  charge  density  and  the  error  begins  another  circuit.  In  the  following  discussions,  the 
exact  value  of/evaluated  at  x,t  is  written  f(x,t) ,  while  the  discrete  approximation  is  written  f' . 


Figure  2  Electrostatic  PIC  error  flow  schematic. 


First,  consider  the  well-known  error  terms  for  the  standard  algorithms.  The  truncation  error 
introduced  by  a  center  difference  of  the  Poisson  Equation  is  well  known  to  be  second  order  in 
space: 

Ax2  w  12  w  360  w  V  ' 

Similarly,  the  truncation  error  in  the  common  leapfrog  difference  for  the  Lorentz  Equations  is 
known  to  be  second  order  in  time  for  the  acceleration: 
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It  is  interesting  to  note  that  the  error  for  x,+At  in  Eq.  16  (found  by  solving  for  x'+A' )  is  higher 
order  than  the  local  error  in  the  central  difference  for  the  velocity,  given  by: 

ev  =  ~ — 7 — “ — ^-x(t  +  At/2)  =  ^-x{3)(t  +  At/2)+-^—  x(S)(t  +  At/2)+o(At6).  17 

At  dt  24  '  1920  V  ’ 

This  is  a  consequence  of  the  cancellation  of  errors  in  the  leapfrog  method,  and  underscores  the 
importance  of  proper  time  centering  of  the  position  and  velocity.  The  error  in  position  for  the 
properly  centered  leap  frog  scheme  is: 

ex  =x,+At  -x(t  +  At)  =  At2ea .  18 


The  charge  density  is  computed  by  interpolating  the  positions  of  the  particles  to  the  mesh.  This 
error  has  been  analyzed  for  a  number  of  interpolation  schemes.  For  nearest  grid  point  weighting, 
the  error  in  charge  density  is  given  by: 

ep,NGp=Pj,NGP -p(x)=^P(2)(x)+-^Pi4){x)+0( Ax6).  19 

For  linear  weighting,  the  error  is  given  by: 

ePMn  =  Pj,lin-p{x)  =  ^-pi2)(x)+^p{4)(x)+o(Ax6).  20 

It  is  especially  illuminating  to  note  that  since  p(2)(x)  =  -£d>(4)(x)  by  differentiating  the  Poisson 
Equation,  the  truncation  errors  in  the  center  difference  of  the  Poisson  Equation  are  exactly 
cancelled  by  the  error  for  linear  weighting  of  charge  density,  to  all  orders  of  accuracy.  Higher 
order  weightings  do  not  have  this  property;  it  is  precisely  cancellation  of  error  that  we  seek  in 
constructing  these  coupled  algorithms,  so  certain  methods  which  are  lower  in  accuracy  for  one 
error  term  often  cancel  other  errors,  resulting  in  superior  accuracy  of  the  coupled  set.  We  shall 
see  below  that  linear  weighting  still  introduces  errors  in  the  interpolation  of  electric  field  to 
particle  location. 

For  quadratic  charge  interpolation,  the  error  is  given  by: 

ep,quad  =  P  j^uud  -  p(x)=  p(2)  {x)+^~  p{4)(x)+  o(Ax6).  21 

o  lO 

The  electric  field  on  the  mesh  must  be  computed  by  finite  differencing  the  potential  on  the  mesh. 
The  difference  is  typically  a  two-cell  center  difference,  in  which  the  error  is  given  by: 

eE  =  °X~A2^X+AC  -  E(x)=  -^-<b(3)(x)+  o(Ax4  ).  22 

The  final  errors  to  consider  in  the  electrostatic  case  are  the  errors  of  interpolation  of  fields  to 
particle  positions.  Although  we  only  consider  the  electric  field  here,  a  similar  exercise  applies  for 
static  magnetic  fields.  The  error  for  the  nearest  grid  point  interpolation  is  given  by: 
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=  /Ax<D(2)(x  +  fAx)-  1  +  3j^  A*V3)(x  +  /Ax)+  o(Ax3 ) 
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where  -l/2</<l/2  is  the  fraction  of  a  cell  width  from  a  mesh  on  which  the  electric  field  is 
known.  Note  that  NGP  weighting  has  a  first  order  error  term. 

For  linear  interpolation,  the  field  interpolation  error  becomes: 

eEi ,  =  (l  -  -  E(x  +  fAx) 

ElM  V  7  ’  2Ax  J  2Ax 

2  24 

=  ^  ~3^~1Ax2g>(3)  (x  +  /Ax)+  o(Ax3 ) 

where  here  0  <  /  <  1 .  Note  that  linear  weighting  is  second  order  accurate. 


For  quadratic  interpolation,  the  error  becomes: 
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^  Ei, quad 


(1  -  2/ ya>. 


2Ax 
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■E(x  +  fAx), 
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=  -7+^/2Ax20(3)  (x  +  /Ax)+  0(Ax3 ) 


where  -  l/2</<l/2.  Note  that  the  quadratic  interpolation  is  second  order  accurate. 


Next,  the  above  analysis  is  extended  to  the  fully  relativistic  electromagnetic  PIC  algorithms.  The 
error  flow  schematic  is  shown  in  Figure  3.  Note  that  the  magnetic  field  is  often  solved  in  parts 
(Faraday  box)  by  advancing  exactly  Vi  to  get  to  a  time  level  t  for  use  in  the  particle  interpolation. 
The  remaining  half  advance  is  performed  just  prior  to  the  Ampere  box.  It  is  useful  to  note  that 
although  several  parts  of  the  diagram  have  multiple  branches,  many  of  the  error  terms  are 
identical  to  the  electrostatic  case.  The  interpolation  of  fields  to  particle  locations  is  unchanged 
from  the  preceding  analysis  (although  the  locations  may  change  by  half  cells  to  accommodate  the 
Yee  mesh  or  an  additional  averaging  of  fields  to  nodes,  shown  above,  may  be  more  efficient). 
The  analysis  for  the  Lorentz  Equations  is  also  still  valid,  with  the  velocity  v  replaced  by  the 
normalized  momentum,  u  =  yv .  An  initial  analysis  of  the  acceleration  error  indicates  the  leading 

term  is  second  order  in  At ,  and  is  also  proportional  to  y6 . 
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Figure  3.  Electromagnetic  PIC  error  flow  schematic. 


In  addition  to  the  divergence  equations  examined  above,  we  must  now  consider  the  Maxwell  curl 
equations.  The  field  components  are  defined  in  the  usual  way,  with  locations  as  shown  in  Figure 
4.  Note  that  the  current  density  is  co-located  with  the  electric  field.  The  error  in  Faraday’s  Law 
can  be  written: 
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f  fields  on  the  standard  Yee  mesh. 
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Similarly,  the  error  in  Ampere’s  Law  can  be  written: 
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The  positions  of  the  analytic  terms,  including  the  derivatives,  in  Eqs.  26-27  are  at  the  Yee  mesh 
locations  as  per  Figure  4.  Both  curl  equations  are  second  order  in  space  and  time. 


The  final  error  term  to  consider  involves  weighting  the  current  to  the  mesh  for  the  source  term  to 
the  Maxwell  Equations.  For  symmetric  current  weightings,  such  as  nearest  grid  point  in  all 
directions,  charge  is  not  conserved  on  the  mesh  assuming  a  symmetric  charge  density  weighting. 
Charge  dipoles  grow  in  time  for  such  a  weighting,  leading  to  catastrophic  failure;  methods  exist 
to  reduce  this  problem8.  The  errors  for  symmetric  current  weightings  are  identical  to  the  charge 
density  weighting  of  the  same  type.  In  charge  conserving  current  weighting  schemes,  the  order  of 
the  interpolation  in  the  direction  of  the  current  component  is  one  order  lower  than  the  transverse 
interpolation.  Here,  we  will  only  analyze  the  error  for  the  NGP-linear  current  weighting: 

<v.  »  -  K  (*  ■ +  A*  / : 2)  = jf»  (x  +  Ax/2)+  (*  +  A*/  2) 

a  2  ,28 

+  —  (37  f’0)  (x  +  Ax  /  2)+  Ay  2/ <2’2)  {x  +  Ax  /  2))+  o(Ax3  ) 

72 


2  .  29 

+  (y  +  Ay  /  2)+  Ax2/;2-8  (y  +  Ay  /  2))+  o(a*s  ) 

/  4 U4 

The  above  error  analysis  demonstrates  that  despite  incomplete  analysis  in  the  literature,  it  is 
possible  to  write  truncation  error  terms  for  each  component  of  the  standard  PIC  method  in  both 
the  electrostatic  and  electromagnetic  limits.  We  believe  this  analysis  is  an  invaluable  contribution 
to  the  existing  body  of  knowledge,  which  provides  a  theoretical  basis  for  the  PIC  algorithms  and 
method  as  a  whole. 
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Porting  XOOPIC  from  XGrafix  to  the  Cross-Platform  Qt  Windowing  Library 

A  large  subset  of  the  existing  XOOPIC  GUI,  which  runs  only  under  Unix,  has  been  ported  to  a 
new  cross-platform  windowing  library,  called  Qt,9  that  runs  equally  well  on  Microsoft  Windows 
95/98/NT  and  all  versions  of  Unix.  Figure  5  shows  a  screen  capture  of  the  new  GUI  running  a 
klystron  simulation  on  a  PC.  The  Run/Stop,  Step  and  Open/Close  Diagnostics  buttons  on  the 
main  control  panel  (seen  in  the  upper  left  of  Figure  5)  have  been  implemented. 


Qt  OOPIC 

Modifications  copyrighted  by  Tech-X  Corporation 
http://www.techrfiorne.corn 


Starting  QtOOPIC  version  1.0 


XOOPIC  (c)  Copyright  1 995-98  The  Regents  of  the  University  of  California 
Plasma  Theory  and  Simulation  Group 
University  of  California  at  Berkeley 
[H  http://ptsg.eecs.berkeley.edu 
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Figure  5.  OOPIC  simulation  of  a  klystron  with  cylindrical  symmetry,  running  under 
Windows  98.  Control  and  diagnostic  windows  are  on  the  left;  top  center  shows  metal 
structures  in  the  z-r  plane  and  electrons  streaming  in  the  z-direction;  bottom  center  shows 
bunching  of  electrons  in  z-pz  phase  space;  messages  from  the  code  appear  in  the  window 
on  the  right. 


The  3-D  plots  of  the  existing  XOOPIC  GUI  have  been  reimplemented  in  the  new  cross-platform 
prototype  using  the  Qt  library's  direct  support  for  OpenGL,10  which  is  now  the  standard  for 
interactive  3-D  computer  graphics.  Figure  6  shows  a  screen  capture  of  the  new  GUI  running  a 
simulation  of  electrons  streaming  into  a  cylindrical  can  and  ionizing  neutral  Argon.  The  two 
lower  right  windows  show  3-D  surface  plots  generated  with  OpenGL.  The  first  three  sliders 
along  the  left  side  of  these  windows  allow  the  user  to  rotate  the  plot  about  the  three  axes,  while 
the  fourth  slider  zooms  the  image  in  and  out.  Axis  labels,  color  effects  and  other  features  have 
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not  yet  been  added.  The  prototype  executable  that  is  available  from  the  Tech-X  web  site  does  not 
include  the  OpenGL  graphics,  because  these  have  not  been  tested  as  thoroughly  as  the  2-D  plots. 


Figure  6.  OOPIC  simulation  of  electrons  streaming  into  a  cylindrical  can  of  Argon, 
including  collisional  ionization.  Top  center  shows  electron  beam  in  the  z-r  plane, 
expanding  radially  due  to  space  charge  forces;  top  right  shows  z-pz  phase  space,  with 
streaming  electrons  along  the  top,  cold  ions  along  the  pz=0  axis,  and  scattered  electrons 
throughout;  center  left  shows  the  number  of  electrons  (top  line)  and  ions  (bottom  line)  in 
the  simulation;  bottom  center  shows  OpenGL  surface  plot  of  the  radial  electric  field 
component  in  the  z-r  plane;  bottom  right  shows  OpenGL  surface  plot  of  the  total  particle 
density 

This  prototype  cross-platform  GUI  has  been  tested  successfully  on  PC's  running  Microsoft 
Windows  95  and  98  and  NT,  as  well  as  on  the  following  Unix  platforms:  a  SunUltra5  running 
Solaris,  a  Pentium  Pro  PC  running  linux,  a  DecAlpha  running  linux  and  an  SGI  running  Irix.  We 
have  established  a  truly  cross-platform  development  environment,  which  allows  us  to  further 
develop  the  GUI  on  any  PC  or  any  Unix  computer,  then  recompile  the  changes  on  all  other 
platforms  with  no  difficulty. 
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